Transient dynamics of molecular devices under step-like pulse bias 
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We report first principles investigation of time-dependent current of molecular devices under 
a step-like pulse. Our results show that although the switch-on time of the molecular device is 
comparable to the transit time, much longer time is needed to reach the steady state. In reaching 
the steady state the current is dominated by resonant states below Fermi level. The contribution 
of each resonant state to the current shows the damped oscillatory behavior with frequency equal 
to the bias of the step-like pulse and decay rate determined by the life time of the corresponding 
resonant state. We found that all the resonant states below Fermi level have to be included for 
accurate results. This indicates that going beyond wideband limit is essential for a quantitative 
analysis of transient dynamics of molecular devices. 
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Anticipating a variety of technological applications, 
molecular scale conductors and devices are the subject 
of increasingly more research in recently years. One of 
the most important issues of molecular electronics con- 
cerns the dynamic response of molecular devices to ex- 
ternal parametersii^i^i^i^i^ii. For ac quantum transport 
in such small devices, atomic details and non-equilibrium 
physics must be taken into account. So, in principle, one 
should use the theory of non-equilibrium Green's function 
(NEGF)^ coupled with the time-dependent density func- 
tional theory (TDDFT)^ to study the time-dependent 
transport of molecular devices. Practically, it is very dif- 
ficult to implement it at present stage due to the huge 
computational cost. One possible way to overcome this 
problem is to use the adiabatic approximation, an ap- 
proach widely used in mesoscopic physics. In this ap- 
proach, one starts from a steady-state Hamiltonian and 
adds the time dependent electric field adiabatically. This 
is a reasonable approximation since most of the time the 
applied electric field is much smaller than the electro- 
static field inside the scattering region. In addition, it 
has been shown numerically® that dc transport properties 
such as I-V curve obtained from the equation of motion 
method coupled with TDDFT agrees with that obtained 
by the method of NEGF coupled with the density func- 
tional theory (DFT) ^°i^^ . Hence, under the adiabatic 
approximation, one could replace TDDFT by DFT and 
use the NEGF+DFT scheme to calculate ac transport 
properties of molecular devices. 

We consider a system that consists of a scattering re- 
gion coupled to two leads with the external time depen- 
dent pulse bias potential Va{t). The time-dependent cur- 
rent for a step-like pulse has been derived exactly going 
beyond the wide-band limit by Maciejko et al^. Since 
the general expression for the current involves triple in- 
tegrations, it is extremely difficult if not impossible to 
calculate the time-dependent current for real systems like 
molecular devices. In this regard, approximation has to 
be made in order to carry out time-dependent simula- 
tions of molecular devices. We note that the simplest 
approximation is the so called wide-band approximation 



where self-energies S""'" are assumed to be independent 
of energy .~ Indeed, if such an approximation is used, i.e., 
S*" = —iT/2, one recovers the expression of transient cur- 
rent first obtained by Wingreen et a,&. However, there 
are two problems when applying this approximation to 
investigate the dynamics of molecular devices. First of 
all, one assumes implicitly that the contribution to the 
transient current is dominated by only one resonant level 
with a constant linewidth function F in the system in such 
an approximation. As we shall show below that this is 
not a good assumption in first principles investigation of 
the dynamics of molecular devices because there could 
be several resonant levels that significantly contribute to 
the transient current in molecular devices. Secondly, in 
the steady state limits at t = and = c» the wide-band 
limit can not reproduce the correct dc I-V curve obtained 
from first principles. By assuming the wide-band limit 
one can get a very different current that depends on the 
choice of F. In this paper, we propose an approximate 
formula of transient current that is suitable for numerical 
calculation for real molecular devices. Our scheme is an 
approximation of the exact solution of Maciejko et al^ 
while keeping essential physics of dynamic systems. Us- 
ing this scheme, we have calculated the transient current 
for several molecular devices. We found that all the reso- 
nant states below Fermi level contribute to the transient 
current. Each resonant state gives a damped oscillatory 
behavior with frequency equal to the bias of pulse and 
decay rate equal to its life time. Because of sharp reso- 
nances, it takes much longer time for the current to relax 
to the equilibrium value. For instance, for a. Al — Ci~ Al 
structure with a transit time of L/vp = 1.3/s the relax- 
ation time is about 50fs. For a CNT-DTB-CNT structure 
with a transit time of Ifs, however, the relaxation time 
can reach several ps due to the resonant state with long 
lifetime. Our results indicated that going beyond wide- 
band limit is crucial for accurate predictions of dynamic 
response of molecular devices. 

From Refilil, the current is expressed as [h = q = 1), 
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(1) 



where 



Mt, e) = A„{t, e)S<'°(e) + ^ ^/^(t, e)S<'"(e)i^^o(t, e)(2) 

where E^''^ and S"'*^ are equiUbrium self-energies and 
Aa{t,e) and Fpa{t,£) have different definitions for up- 
ward and downward pulses (see Ref 4 for details). In the 
absence of ac bias, Aa{t, e) is just the Fourier transform 
of retarded Green's function. As discussed in RefH that 
the first term in Eq.Q corresponds to the current flow- 
ing into the central scattering region from lead a while 
the second term corresponds to the current flowing out 
from the central region into lead a. From Eq.([T]) we see 
that in order to calculate the transient current for a pulse 
bias we need to include the states with energy from —oo 
to the Fermi energy. This is very different from dc case 
where only the states with energy in the range vl — vr 
about Fermi level contribute. Physically, this can be un- 
derstood as follows. For ac transport with a sinusoidal 
bias cos{Lot), the photon assisted tunneUing is significant 
only for the first a few sidebands^. The step- like pulse 
can be expanded in terms of sinusoidal bias with contin- 
uous distribution of frequencies and each sinusoidal bias 
generates a photon sideband that facilitates the photon 
assisted tunnelling. Hence we expect that all the resonant 
states below Fermi level should be included and carefully 
examined in the calculation of transient current. Note 
that Eq.((T]) and ([2]) are exact expressions with Aa{t,e) 
and Ff}a(t, e) given in Ref 4. Our approximation is made 
on Aa{t, e) and Fpa{t, e). For the upward pulse, Aa{t, e) 
and Fj3a{t, e) are given by the following ansatz. 



A^(i,e)=A5'Ji,e)+AL(t,e) (3) 



with 



and 

e)]t = E50(e)A?^ + E^f{e - + vp)A-^p (5) 

where 



(6) 



Gl{E,e) = ll[E-H-U-Y, S^'"(e + v^- vp)] (7) 

with U(,q and U are, respectively, the equilibrium 
Coulomb potential and dc Coulomb potential at bias 



Vl ^ vji. As will be illustrated in the examples given be- 
low this ansatz can be easily implemented to calculate the 
transient current for real molecular devices. Importantly, 
the results obtained from the ansatz captured essential 
physics of molecular devices. We wish to emphasize that 
our ansatz goes beyond the wide-band limit. It agrees 
with the expression of time-dependent current obtained 
by Wingreen et al in the wide-band limit^'^ and produces 
correct limits at t = and t —>■ oo. 

Note that Gq{E, e) and G^iE, e) are different from the 
usual definition of Green's functions, they allow us to 
perform contour integration over energy E in Eq. ^ and 
([5]) by closing a contour with an infinite radius semicir- 
cle at lower half plane. For a constant e, we have the 
following eigen equations 



(8) 



Expanding G5(i?, e) in terms of its eigen functions \ipn) 
and |(/>^j), we havei^ 



GS(i?,e)=5]|V^°)(0°|/(i?-6°+»O+) 



(9) 



With similar expression for G^{E, e), Eq.(|4]) can be writ- 
ten asi^ 



1 _ pj(e-e,i+t)Q)t 

n 

Now we show that our formalism gives the correct lim- 
its. At t = we have A^(i,e) = G5(e) and F^Jt,e) = 
Gg(e)I]^'°(e) with G5(e) the equihbrium Green's func- 
tion. This shows that the current from Eq.(l2]) is zero. 
Since all the poles e° and e„ in Eq. pU]) are on the lower 
half plane, at i ^ oo we have = G^ {e + Va) and 
Fpc. = + «/3)S°'°(e + vp - v^) where G'-(e) is the 
Green's function with dc bias Va ai t oo. Substi- 
tuting expressions of A^ and F^^ into Eq.([T]), it gives 
the same dc current at the bias — v^. So far, we 
have discussed the ac conduction current Ja(t) under 
pulse-like bias. The displacement current due to 
the charge pileup dQ/dt inside the scattering region can 
be included using the method of current partitio n^^i^^ : 
— ~{Jl + Jr)/2, so that the the total current is given 
by h = {Jl - Jr)/2& 

With the formalism established, we now proceed to 
calculate the dynamic response of molecular devices. We 
have used the first principle quantum transport package 
MatDcal.^" To calculate the transient current for step- 
like pulse, we need to go through the following steps: 
(1). calculate two potential landscapes using NEGF- 
DFT package: the equilibrium potential Ueq at i = 
and the dc potential U at t — oo. (2). With Ueq and 
U obtained, one solves eigenvalue problem using Eq.® 
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FIG. 1: (color online) Transient current of AI — C4 — Al struc- 
ture at different bias vl = 0.001, 0.003, 0.01, 0.03a«. The solid 
line shows the transient current and the dotted line is the dc 
current Idc at bias vl — «h = 2vl- Inset of Fig. la: schematic 
plot of the Al — C4 — Al structure. 



FIG. 2: (color online) Transient current of CNT-DTB-CNT 
structure at different bias vl = 0.001, 0.003, 0.01, 0.02a«. The 
solid line shows the transient current and the dotted line is the 
dc current at bias vl — vr = 2vl- Inset of Fig. 2b: schematic 
plot of the CNT-Di-thiol benzene-CNT structure. 



and its counterpart for U, then find A^^ and from 
Eg. prHl . and finally A'^ and Fg^ can be calculated from 
Eq.dS]) and p. 

Inset of Fig. la shows the structure of AI—C4—AI where 
Al leads are along (100) direction. The nearest distance 
between Al leads and the carbon chain is 3.781au and 
the distance of C-C bond is 2.5au.{lau — 0.529A). Fig.l 
shows the total transient currents lL{t) of the AI — Ca~AI 
structure with various voltages wl(^) = —VR{t). Follow- 
ing observations are in order. First of all, for all bias 
voltages the transient currents reached the correct limits 
at i = and t = 00. Secondly, we see that once step- like 
voltage is turned on in the lead, currents oscillate rapidly 
with large amplitude in the first a few fs and then grad- 
ually approach to the steady-state values {Idc shown in 
the figure). In the first 10 to 30 fs, the current is much 
larger than that of the steady state value which agrees 
with the results obtained using first principle calculation 
with wide-band limit^. For Fig. la, the relaxation time 
(time to reach to steady state) is roughly 150 fs and for 
Fig.lb-ld the relaxation time is about 50fs. In addition, 
the switch-on time (the time to reach the maximum cur- 
rent) for the Al — C4 — Al structure does not depend on 
the bias voltages. The typical switch-on time is about 
2fs for applied bias voltage vl ranging from O.OOlau to 
O.Olau {lau = 27. 2^). For Al leads, the Fermi veloc- 
ity is about 2 x 10^ m/s which corresponds to a transit 
time of 1.3 fs for the Al — C4 — Al structure whose size is 
about L ^ A7 au. Thirdly, we observe that the dc limit 
Idc at Vl — O.Olau is larger than that at — 0.03au. 
This is due to the appearance of the negative differential 
resistance at about vl = 0.0075au. Finally, there are 
several timescales characterizing the dynamic response 
of the molecular device. This can be seen clearly from 
Fig. Id that after 10 fs, the system shows a damped oscil- 
lation similar to the charging process of a classical RLC 



circuit. We will discuss this kind of oscillation in detail 
in the second example. 

As a second example, we study the transient current 
for di-thiol-benzene molecule (DTB) in contact with two 
(3,3) carbon nanotube (CNT) leads (see inset of Fig. 2b). 
The structure is relaxed with the distance between the 
S atom and the nearest C atom equal to 2.73aM and 
the bond length of C-C being 3.61aM. Fig. 2 shows the 
transient current for different upward pulse biases. We 
see that for small bias vl = O.OOlau, the current drops 
quickly in first 50 fs and then oscillates with much slower 
decay rate. It is found that the oscillatory part of the 
transient current is dominated by cos(2vLt) which re- 
mains valid for the transient current at other biases vl 
shown in Fig. 2b to Fig. 2d. For instance, this gives the 
distance between adjacent peaks tq = 76/s in Fig. 2a 
when VL = O.OOlau. With lau = 0.00242/s, we obtain 
To — tt/vl- Different from the Al — C4 — Al structure, it 
takes much longer time for the system to reach the equi- 
librium current /l(oo) = 1.5/i^ (after SOOOfs the current 
is about 1.7fiA). From Fig.2a-2d, we conclude that the 
relaxation time is several ps. 

Physically, this can be understood from the transmis- 
sion coefficient T{E). Fig. 3 depicts T{E) vs energy rang- 
ing from the transmission threshold to Fermi energy. We 
have scanned 100,000 energy points in order to resolve 
sharp resonant peaks labelled in Fig. 3. In our calcula- 
tion, 100 energy points were used for each sharp resonant 
peak (total 3000 energy points used) to converge the inte- 
gration over e, i.e., / dtTT[Ja{t, e)] in Eq.(l). Since these 
sharp resonant peaks correspond to resonant states with 
large lifetimes, the incoming electron can dwell for a long 
time at these resonant states and hence the correspond- 
ing current decays much slower than the other states. 
If we focus on a particular resonant state with resonant 
energy eg (below Fermi level) and half- width Fq, then 
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FIG. 3: (color online) The transmission coefficient for the 
CNT-DTB-CNT structure. The insets are T{E) vs energy at 
peak 7 and its corresponding current. 

Eq.p gives ~ exp[i(e - eo + Va)t - {To/2)t]i2.. As- 
suming that the sharp resonant state gives major contri- 
bution to the current (the wideband approximation), wc 
have A" ~ exp{ivat — (ro/2)i). Therefore the first term 
in Eq.(l2|) exhibits an oscillatory part exp{ivat — (ro/2)i) 
while the second term behaves like ex-p{2iVat — rot). It is 
the interplay between these two terms that gives rise to 
the transient current. For instance, for Fig. 2 the second 
term exp{2ivat — Tot) dominates while for Fig. Id the first 
term gives the most contribution. 

Indeed, our numerical result confirms this analysis. It 
shows that these resonant peaks give major contributions 
to the transient current for t > 50/s. In addition, we find 
that there is an one to one correspondence between the 
resonant peak at e and the corresponding Tr[j7Q(i, e)]: 
Tr[^7Q(t,e)] exhibits a huge peak whenever e is near the 
resonance. This correspondence is important because it 
indicates that our ansatz has kept essential physics aris- 
ing from the above analysis. Furthermore, our result 
shows that the transient current due to each resonant 



peak has the same characteristics frequencies or 2vl- 

Let's examine the contribution of each resonant state 
to the oscillatory part of the transient current at vl — 
O.OOlau (Fig. 2a). Among these resonant peaks in Fig. 3, 
the most contribution comes from the peak number 7 
with half-width r7 = 2.7 x lO^^au which corresponds to 
a decay time ry = 1240/s from the expression exp(— r7i). 
In the left inset of Fig. 3, we plot the current obtained by 
integrating Tr[j7Q(e)] over the neighborhood of peak 7 
(see right inset of Fig. 3). It shows that the decay time is 
indeed characterized by tj. Comparing Fig. 2 and Fig. 3, 
we see that the contribution from the peak 7 to the total 
current is about 15% for t < 50/s while for t > 50/s the 
contribution is 50%. The next dominant contribution is 
due to the peaks numbered 5,10, and 12 whose contribu- 
tions are one order of magnitude smaller. This indicates 
that one has to include all the resonant peaks for accu- 
rate results. Since different resonant peak corresponds 
to a different half- width F, one can not choose just one 
r to characterize the system. We have also calculated 
the transient current for the structure of Al — Cgo ~ Al 
and our results show that the long time behavior is dom- 
inated by two resonant peaks with different F and shows 
beat pattern with relaxation time about 800fs. 

In summary, we have carried out first principles inves- 
tigation of time response of molecular devices. We found 
that the resonant states below Fermi level are crucial for 
time-dependent quantitative analysis. Our results indi- 
cated that the long time behavior of transient current 
is dominated by resonant states and the individual res- 
onant state gives the damped oscillatory behavior with 
frequency equal to the bias of pulse and decay rate equal 
to the life time of the corresponding resonant state. Our 
results indicated that one has to go beyond the wide-band 
limit for quantitative calculations of dynamic response of 
molecular devices. 
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In ab initio calculation, if atomic orbits are chosen, their 
overlap matrice S have to be included in all Green's func- 
tions by replacing E with ES. In this case Eq.® has to 
be divided by a factor < (;/)5^|5'|i/)5J >. 

In order to use this formalism (Eq. (|10|l ') correctly, one must 
be careful in the self-energy calculation. In the numerical 
calculation of self-energy, one usually adds a small imag- 
inary part to the real energy E —> E + i-q in order to re- 
solve the retarded or advanced self-energies. Unfortunately, 
this parameter r; will introduce many spurious poles in the 
lower half plane with imarginary part of the poles less than 
rj. This has no effect on ac current if it is calculated directly. 
However, if the theorem of residue is used such as Eq. pop 
the transient behavior will be dominated by the spurious 
poles. To eliminate this effect, one has to calculate the 
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